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I  Introduction: 


The  implementation  of  inverse  planning  technique  is  essential  to  the  realization  of  the  potential  of 
intensity  modulated  radiation  therapy  (IMRT).  The  most  popular  inverse  planning  techniques  currently 
used  fall  into  two  distinct  categories:  1)  beamlet-based  inverse  planning  and  2)  aperture-based  inverse 
planning.  In  beamlet-based  planning,  each  field  (a  gantry  angle  -  energy  combination)  is  divided  into  a 
matrix  of  beamlets  whose  weights  are  optimized  independently.  The  optimization  therefore  yields  an 
intensity  map  which  is  then  is  translated  into  a  sequence  of  leaf  aperture  shapes  before  being  delivered. 
There  are  several  problems  associated  with  beamlet-based  inverse  planning.  First,  optimized  intensity 
maps  are  often  converted  into  a  series  of  discrete  levels  in  an  attempt  to  make  the  sequencing  easier. 
However  this  step  also  introduces  some  quantization  errors  and  therefore  results  in  loss  of  treatment  plan 
quality.  Second,  leaf  sequencing  is  constrained  by  hardware-related  factors  and  therefore  often  requires  a 
large  number  of  complex  field  shapes  to  deliver  a  given  intensity  map,  thereby  decreasing  the  overall 
efficiency.  Third,  since  the  leaf  sequencing  step  is  excluded  from  the  intensity-map  optimization  therefore 
all  the  delivery-related  effects  such  as  leakage,  the  tongue-and-groove  design  and  head  scatter  are  not 
taken  into  account  when  choosing  an  intensity  map.  Efforts  to  include  these  effects  have  also  introduced 
significant  additional  complexity.  Aperture-based  optimization  is  designed  to  reduce  the  complexity  of 
IMRT  treatment  plans  and  is  flexible  enough  to  easily  include  delivery-related  effects.  In  particular, 
aperture  optimization  requires  no  leaf  sequencing  and  apertures  are  guaranteed  to  satisfy  delivery 
constraints.  It  is  therefore  a  promising  alternative  to  beamlet-based  IMRT.  In  this  work,  we  consider  a 
method  of  contour-based  treatment  planning  in  which,  unlike  in  Direct  Aperture  Optimization  (DAO),  the 
aperture  shapes  are  initially  chosen  to  be  the  Beam  Eye’s  View  (BEV)  of  the  target  and  the  critical 
structures,  and  therefore  are  given  as  constants  in  the  optimization  process.  The  obvious  advantage  of  this 
method  with  regard  to  DAO  is  that  the  only  parameters  that  need  be  optimized  are  the  fields’  weights,  and 
the  plan  is  not  optimized  over  all  the  possible  aperture  shapes.  Our  hope  is  naturally  that  optimal  aperture 
shapes  are  somewhat  close  to  the  BEV  and  that  making  this  choice  will  save  a  significant  amount  of 
computation  time. 

II  Research  and  Accomplishments 

We  now  turn  to  the  description  of  the  optimization  algorithm.  The  point  of  the  following  optimization 
scheme  is  to  find  the  beamlet  weights  that  will  deliver  a  dose  conformation  that  is  as  close  to  the 
prescription  as  possible.  A  natural  way  is  to  use  an  iterative  approach  in  order  to  find  the  optimal  plan:  an 
initial  plan  is  evaluated,  adjustments  are  made,  the  plan  is  reevaluated,  further  adjustments  are  made,  and 
so  on  until  some  convergence  criterion  is  met.  In  order  to  evaluate  a  treatment  plan,  we  need  to  be  able  to 
assess  ‘how  far’  it  is  from  the  prescription.  This  evaluation  can  be  made  by  condensing  the  quality  of  a 
plan  into  a  single  value.  This  value  is  referred  to  as  the  objective  function  and  naturally  depends  on  the 
choice  of  criteria  that  was  made  to  compute  it.  The  objective  function  is  therefore  a  mathematical  function 
that  takes  as  its  input  the  dose  distribution  of  both  the  evaluated  plan  and  the  prescription.  Once  we  have 
determined  the  criteria  that  define  the  objective  function,  we  can  easily  compare  two  possible  treatment 
plans  by  comparing  the  values  of  their  respective  objective  functions. 

In  the  following  approach,  the  objective  function  is  chosen  to  be  roughly  the  sum  of  the  squared 
differences  between  the  prescription  and  the  delivered  doses  at  different  measurement  points  in  different 
structures.  We  naturally  want  to  achieve  a  plan  whose  objective  function  is  as  small  as  possible  (since  a 
small  objective  function  indicates  that  the  delivered  doses  are  close  to  the  prescription  ones).  In  order  to 
derive  the  expression  of  the  objective  function,  we  first  need  to  introduce  some  notations. 
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Figure  1 :  Comparison  of  flowcharts  for  beamlet-based  and  aperture-based  IMRT.  The  segment  weigths 
optimization  is  not  detailed  here  for  clarity. 

2.1  Notations  for  Plan  Size  Parameters,  Doses  and  Constraints 

We  begin  with  a  CT  dataset  containing  structures  of  interest  that  are  manually  contoured.  There  are  A 
structures  (organs).  Index  a=0  refers  to  the  target  structures  while  other  structures  are  indexed  from  1  to 
A-l.  Each  structure  a  contains  Na  measurement  points,  for  a=0,K,^-l.  We  also  assume  that  organs  do  not 
overlap,  which  means  that  every  measurement  point  belongs  to  one  and  only  one  structure. 

As  we  discussed  earlier,  beam  geometries  and  energies  are  specified  before  the  optimization.  Each  (beam 
angle/energy)  combination  is  called  a  field;  a  port  refers  to  a  single  beam  geometry.  Just  like  in  beamlet- 
based  optimization,  each  field  is  divided  into  beamlets,  and  we  denote  B  the  total  number  of  beamlets.  All 
the  beamlets  belong  to  a  total  number  of  G  fields  (every  beamlet  only  belongs  to  one  single  fields).  bg  is 
the  number  of  beamlets  in  field  g,  and  we  thus  have  that  B  =  . 

f'  1 

For  each  beamlet  b  (where  b  e  [l;fi])  and  each  structure  a,  we  define  the  dose  kernel  as  the  array  of 
doses  received  by  each  point  n  in  structure  a  when  only  beamlet  b  is  given  weight  1  and  all  the  other 
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beamlets  are  given  zero  weights.  Since  structure  a  contains  Na  points,  we  have  d{ba)  =  {af£?,K  ,«/$  }•  Using 

accurate  values  for  the  beamlet  kernels  is  critical  to  the  final  treatment  plan  quality,  and  this  is  why  dose 
kernels  are  computed  via  Monte-Carlo  simulation  in  this  work. 

The  weight  given  to  beamlet  b  is  denoted  wb,  and  we  are  trying  to  determine  the  set  of  weights  that  best 
matches  the  prescription.  Once  the  dose  kernels  are  computed  (and  they  are  only  computed  once  prior  to 
the  optimization),  the  doses  received  at  each  point  in  each  structure  at  any  time  depend  only  on  the 
beamlet  weights  at  that  time.  This  is  why  the  vector  w  =  {w,,K  ,we}  that  contains  all  the  beamlet  weights 
is  referred  to  as  the  state- vector. 

Indeed,  knowing  the  dose  kernels  and  the  state  vector,  we  can  easily  compute  the  accumulated  dose  at 
every  measurement  point  (or  voxel).  We  denote  D f 1  as  the  total  accumulated  dose  at  point  n  in  organ  a. 

In  particular,  we  have  that  n<°>  =  V  w  du,) ,  which  is  a  just  a  mathematical  way  to  express  that  the  total 

r  n  vb**b,n 

b= 1 

accumulated  dose  at  point  n  in  organ  a  is  simply  the  sum  of  the  doses  delivered  by  each  individual 
beamlet. 

We  now  need  to  specify  our  conventions  regarding  prescription  doses,  which  are  defined  as  follows.  A 
prescribed  dose  is  specified  for  the  target,  and  then  dose-volume  constraints  are  defined  for  both  the  target 
and  critical  structures.  We  denote  d<°>  as  the  prescribed  dose  for  the  target.  Besides  this  prescribed  dose, 

there  are  two  additional  dose-volume  constraints  that  apply  to  the  target.  The  upper  dose-volume 
constraint  is  of  the  form:  “no  more  than  v,(0)%  of  organ  0  (i.e.  the  target)  should  receive  a  dose  greater 

than  Z)l(0)  ”,  whereas  the  lower  dose- volume  constraint  is  of  the  form  “no  more  than  v,(0)%of  organ  0 
should  receive  a  dose  less  than  D,(0)  ”. 

For  each  other  structure  a,  there  are  Ca  dose-volume  constraints,  each  one  of  them  of  the  form  “no  more 
than  v*a)%  of  organ  a  should  receive  a  dose  greater  than  £)'a)”.  Eventually  each  constraint  c  that  apply 
to  organ  a  (including  the  target)  is  given  a  relative  weight  factor  a(ca) . 

2.2  Objective  Function  and  Gradient  for  Beamlet-Based  IMRT 

In  these  conditions,  the  objective  function  is  the  sum  of  term  that  represents  the  target  prescription  dose, 
plus  other  terms  that  represent  violations  to  all  the  constraints  that  apply: 

a- 0  c=l 

The  term  for  the  target  prescription  dose  has  the  following  quadratic  form: 

It  can  be  shown  (see  Appendix  for  details)  that  the  gradient  of  the  objective  function  has  the  following 
form: 


dF 

dwb 


= -2  ■  -  ^<$>1 

a- 1  ™  a  c,a 


(4) 
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Once  we  have  computed  a  gradient,  we  can  move  along  its  direction  (departing  from  the  initial  state)  and 
find  the  step  size  that  minimizes  the  objective  function.  This  step  is  referred  to  as  the  line-minimization, 
and  comes  down  to  finding  step  size  t  such  that  F(w+/VF)  is  minimal.  We  first  turn  to  the  adaptation  of 
this  method  to  aperture-based  optimization  before  tackling  the  additional  difficulties  that  pertain  to  the 
line  minimization  step. 

2.3  Objective  Function  and  Gradient  for  Aperture-Based  IMRT 


In  aperture-based  optimization,  we  want  to  express  this  function  in  terms  of  the  weights  of  each  field.  In 
other  words,  we  want  to  set  the  weights  of  all  the  beamlets,  within  each  field,  equal.  We  must  also  make 
sure  that  this  simplification  remains  throughout  the  optimization  process,  and  thus,  that  the  gradient 
direction  along  which  we  evolve  imposes  the  same  changes  in  weights  for  all  the  beamlets  that  belong  to 
the  same  field. 


In  order  to  achieve  this,  we  need  to  modify  the  gradient  computation.  We  thus  can  consider  a 
modification  of  expression  (3)  in  which  we  incorporate  the  fact  that  all  beamlets  from  the  same  field  have 
the  same  weights.  This  yields  the  following  expression  of  the  objective  function: 


f = oi‘{  1<:. 


'  N 

a- 1  1  ’  a  c,a  w=l 
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where  G  is  the  number  of  ports,  Bg  is  the  number  of  beamlets  in  port  g. 


(5) 


If  we  duplicate  the  proof  through  which  we  derived  (4)  from  (3),  we  get  a  modified  version  of  expression 
(4)  that  incorporates  the  fact  that  the  objective  function  must  only  be  differentiated  with  respect  to  all  the 
beamlet  weights  that  belong  to  the  same  field.  This  naturally  dramatically  reduces  the  number  of 
differentiations,  thus  limiting  the  dimension  of  the  space  over  which  we  have  to  optimize.  The  new 
gradient  expression  is  now: 


We  must  now  take  into  account  the  effect  of  the  non-negativity  constraints,  which  already  existed,  in  the 
original  process.  These  constraints  simply  account  for  the  fact  that  beamlet  weights  must  remain  non¬ 
negative  since  a  negative  weight  has  no  physical  meaning,  even  though  we  can  conceptually  imagine  that 
‘removing’  dose  from  the  patient  would  reduce  the  objective  function.  Negative  coefficients  can  occur 
for  instance  if  we  depart  from  a  valid  initial  state,  i.e.  such  that  wb  £  0,  for  all  b,  and  then  perform  a  line 
minimization  along  the  gradient  direction,  which  yields  a  step  size  of  t.  If  there  is  a  b  such  that  thb  <  , 

then  the  new  weight  for  beamlet  b  will  become  negative.  In  order  to  avoid  such  aberration,  we  must  find  a 
t  that  is  such  that  thb  >  -wb  for  all  beamlets,  and  this  can  be  done  by  performing  a  loop  through  all  the 

beamlets  and  scaling  t  down  to  -wb/hb  whenever  thb  <  -wb .  By  scaling  the  initial  t  down  to  the  t’  yielded 

by  this  loop,  we  are  actually  guaranteed  (under  technical  yet  general  conditions)  that  the  new  step  will 
result  in  some  reduction  of  the  objective  function  since  we  are  still  moving  in  the  direction  of  the 
minimum.  This  loop  is  made  clear  in  the  following  pseudo-code: 


t’=t; 

for  b  =  1  to  B, 

w’b  =wb  +thb; 
if  w’b  <  0  then 
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t’=  -  wb  /  hb 


end  if 

_ end  for _ _ _ 

Alg.  1 :  Scaling  step  in  order  to  avoid  negative  weights. 

However,  several  problems  arise  from  such  scaling.  If  there  is  a  beamlet  b  for  which  =  0  and 

thb  <  0 ,  then  t  will  be  scaled  down  to  zero.  One  way  to  solve  this  problem  is  to  modify  the  search 

direction  before  looping  through  beamlets,  and  in  particular  to  set  to  zero  the  negative  hb  that  correspond 
to  zero-weighted  beamlets  since  we  know  that  these  beamlets  will  scale  the  step  size  down  to  zero  if 
nothing  is  done  preventively.  However,  by  doing  this  we  modify  the  search  and  we  therefore  have  no 
other  guaranty  that  the  step  size  will  be  negative.  We  therefore  need  to  come  up  with  two  search 
directions  (one  for  positive  steps  t,  one  for  negative  steps  t)  along  which  the  step  size  can  be  scaled  down. 
This  method  is  detailed  in  algorithm  2. 

Eventually,  once  we  have  the  two  search  directions,  we  just  need  to  find  t  that  minimizes  the  following 
function:  q( a  _  frif  +  h  )  if  t<.0 ,  an(j  nutnerica]  techniques  to  bracket  the  minimum  are  employed  but 

W  +  A+)  Vt> 0 

not  discussed  here. 


for  b  =  1  to  B, 

if  (wb  =  0  and  hb  >  0)  then 

K  =  0 

K-K 

else  if  (wb  =  0  and  hb  <  0) 

o 

else 

K  =  h 
K  =  h 

end  if 

_ end  for _ _ _ 

Alg.  1 :  Modification  of  the  search  direction  in  order  to  avoid  zero-length  steps 

The  optimization  algorithm  for  the  simplified  objective  function  was  only  tested  on  a  simple  plan  of  one 
single  port  made  20*12  beamlets.  As  one  would  expect,  the  number  of  steps  in  the  optimization  process 
was  dramatically  reduced,  from  1500  down  to  1  single  step.  This  makes  sense  since,  for  a  single  port,  the 
aperture-based  optimization  is  reduced  to  a  single-variable  function  optimization.  No  further  testing  of  the 
algorithm  could  be  performed  before  expiration  of  the  grant. 
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